# file ratdef.r  ssw 9/27/83     modified by jwt 8/21/84
# ratfor macro definitions for inclusion in ratfor programs
#
# file ratdef.r  ssw 9/27/83
DEFINE(CHAR,character)
DEFINE(DIG0,'0')
DEFINE(MINUS,'-')
define (NEWLINE,char(10))
define (BLANK,' ')
define (TAB,char(9))
define (EOS,char(0))
define (EOF, -1)
define (MPTEST,#)		# to omit tests
define (MPBASE,1073741824)
define (MPBITS,30)		#  2**MPBITS = MPBASE
define (MPDIM,120)
define (MPCDIM,600)
define (BASE,23)		# base for sprp test
define (TRUE,1)
define (FALSE,0)
define (STDIN,5)
define (STDOUT,6)		
define (STDERR,0)		# output:  error messages
# file: mp1.r
#include 'primit.r'
# file primit.r   ssw 9/27/83
# modified 6/16/88 by ssw to allow output to any file in putdec
# now the primitive subroutines:
#
integer function getlin(line,ch)  # read a line from channel ch
character line(*)
integer ch,i,j,fgetc
i=1
repeat {
	j=fgetc(ch,line(i))
	if (j != 0) { # EOF (or i/o error)
		line(i)=EOS
		getlin=j
		return
		}
	if (line(i) == NEWLINE) {
		line(i)=EOS
		getlin=0
		return
		}
	i=i+1
	}
end

subroutine putlin(line,ch) # output line on channel ch
integer ch,i,j,fputc
character line(*)
for (i=1; line(i) != EOS; i=i+1)
	j=fputc(ch,line(i))
return
end

subroutine error(str)
character str(*)
call putc(NEWLINE) # flush stream output
call putlin(str,STDOUT)
stop
end

subroutine putdec(cha,n,w)
character ch(79)
integer itoc,i,j,n,nd,w,cha,fputc
nd=itoc(n,ch,79)
for (i=nd+1; i<=w; i=i+1)
     j=fputc(cha,BLANK)
for (i=1; i<=nd; i=i+1)
     j=fputc(cha,ch(i))
return
end

integer function itoc(int,str,size) # convert integer int to char string in str
integer iabs, mod, d,i,int,intval,j,size  # copied from software tools
character k, str(size), digits(10)              # by Kernighan and Plauger
data digits /'0','1','2','3','4','5','6','7','8','9'/

intval=iabs(int) # get absolute value
str(1)=EOS
i=1
repeat { # generate digits
     i=i+1
     d=mod(intval,10)
     str(i)=digits(d+1)
     intval=intval/10
     } until (intval==0|i>=size)
if (int<0 & i<size) {  # put in the sign
     i=i+1
     str(i)=MINUS
     }
itoc=i-1
for (j=1; j<i; j=j+1) {  # then reverse
     k=str(i)
     str(i)=str(j)
     str(j)=k
     i=i-1
     }
return
end

integer function ctoi(in,i) # convert string at in(i) to integer, increment i
CHAR in(1)
CHAR*10 digits
integer index,d,i
data digits/'0123456789'/
while (in(i)==BLANK)
     i=i+1  # find the first non-blank
for (ctoi=0; in(i)!=EOS; i=i+1) {
     d=index(digits,in(i))
     if (d==0)  # if non-digit
          break
     ctoi=10*ctoi + d - 1
     }
return
end

# MP subroutines MPCOPY, ADD, addl, SUB, MUL1, DIV1, mpchk, MUL, DIV
# UPPER CASE subroutines are for general use.
# lower case subroutines are mostly for private internal use only.

# MPCOPY:  copy A into B     ssw 2/26/83
subroutine mpcopy(a,b)
integer a(1), b(1), i, d
d=a(1)    
MPTEST call mpchk(100,a)
do i=1,d
b(i)=a(i)
return
end

# ADD:  C = A + B    ssw 2/26/83
subroutine add(a,b,c)
integer a(1), b(1), c(1), carry, i, t, ba
data ba/MPBASE/
MPTEST call mpchk(200,a)
MPTEST call mpchk(201,b)
if (a(1)==b(1)) { carry=0
     for (i=2; i<=a(1); i=i+1) {
          t = a(i) + b(i) + carry
          if (t<ba) {carry=0; c(i)=t}  # no carry
          else {carry=1; c(i)=t-ba}    # carry 1
          }
     if (carry==1) {c(1)=a(1)+1; i=c(1); c(i)=1}
     else {c(1)=a(1)}  # set length of c in c(1)
     }
else if (a(1)<b(1)) {call addl(a,b,c)}
else {call addl(b,a,c)}  # b(1) < a(1) in the third case
MPTEST call mpchk(202,c)
return
end

# addl: C = A + B; assume A(1) < B(1);
# for use by ADD only.   ssw 2/26/83
subroutine addl(a,b,c)
integer a(1), b(1), c(1), carry, i, t, ba
data ba/MPBASE/
carry=0
for (i=2; i<=a(1); i=i+1) {
     t = a(i) + b(i) + carry
     if (t<ba) {carry=0; c(i)=t}  # no carry
     else {carry=1; c(i)=t-ba}    # carry 1
     }
# we assume that a(1) < b(1)
# omit next initialization because i just continues
for ( ; i<=b(1); i=i+1) {
     t = b(i) + carry
     if (t<ba) {carry=0; c(i)=t}  # no carry
     else {carry=1; c(i)=t-ba}    # carry 1
     }
if (carry==1) {c(1)=b(1)+1; i=c(1); c(i)=1}
else {c(1)=b(1)}  # set length of c in c(1)
return
end

# SUB:  C = A - B;  assume A >= B    ssw 2/26/83
subroutine sub(a,b,c)
integer a(1), b(1), c(1), borrow, i, t, ba
CHAR*80 smsg
data ba/MPBASE/
smsg = 'negative diff in sub'//char(0)
MPTEST call mpchk(300,a)
MPTEST call mpchk(301,b)
if (a(1) < b(1)) call error(smsg)
borrow=0
for (i=2; i<=b(1); i=i+1) {
     t = a(i) - b(i) - borrow
     if (t<0) {borrow=1; c(i)=t+ba}  # borrow 1
     else {borrow=0; c(i)=t}         # no borrow
     }
# omit next initialization because i just continues
for ( ; i<=a(1); i=i+1) {
     t = a(i) - borrow
     if (t<0) {borrow=1; c(i)=t+ba}  # borrow 1
     else {borrow=0; c(i)=t}         # no borrow
     }
if (borrow==1) {call error(smsg)}
else {c(1)=a(1); i=c(1)
     while (i>2 & c(i)==0) {c(1)=c(1)-1; i=c(1)}
     }   #  condense dimension of c
MPTEST call mpchk(302,c)
return
end

# MUL1:  C = A*B;  B single precision     ssw 2/26/83
subroutine mul1(a,b,c)
integer a(1), b, c(1), carry, d, i, x(2)
CHAR*80 negb
negb = ' negative b in mul1'//char(0)
MPTEST call mpchk(400,a)
if (b>0) goto 10
if (b<0) {call putdec(STDOUT,401,6);call putdec(STDOUT,b,10)
     call putlin(negb,STDOUT); call putc(NEWLINE)}
c(1)=2
c(2)=0
return
10 carry=0
d=a(1)
do i=2,d {
     call muladd(a(i),b,carry,x)  # x = a(i)*b + carry
     carry = x(2)     # carry high order
     c(i) = x(1)      # store low order as next digit
     }
if (carry==0) {c(1)=d}
else {c(1)=d+1; i=c(1); c(i)=carry}
MPTEST call mpchk(402,c)
return
end

# DIV1:  remainder = DIV1(a,b,quotient);  divide  a  by  b
# b single precision      ssw 2/26/83
integer function div1(a,b,q)
integer a(1), b, q(1), i, r, y(2)
CHAR*80 negb
negb = '  b <= 0 in div1'//char(0)
if (b>0) goto 10
call putdec(STDOUT,501,6)
call putdec(STDOUT,b,10)
call putlin(negb,STDOUT)
call putc(NEWLINE)
stop
10 r=0
MPTEST call mpchk(500,a)
for (i=a(1); i>1; i=i-1) {
#     x(2)=r
#     x(1)=a(i)
#     r=divrem(x,b,q(i))  # divide two-digit by one-digit
      call divrem(r,a(i),b,y)
      r=y(2)
      q(i)=y(1)
     }
div1=r
q(1)=a(1)
i=q(1)
while (i>2 & q(i)==0) {q(1)=q(1)-1; i=q(1)}
return
end

# mpchk: check dimension and digits; internal use only   ssw 2/26/83
subroutine mpchk(code,a)  # code identifies caller; we check  a.
integer code, a(1), d, i
CHAR*80 s1,s2,s3,s4,s5
s5 = ' mpchk'//char(0)
s1 = ' dimension < 2 '//char(0)
s2 = ' dimension > MPDIM '//char(0)
s3 = ' negative digit '//char(0)
s4 = ' digit >= MPBASE '//char(0)
d=a(1)
call putlin(s5,STDOUT)
call putdec(STDOUT,code,5)
do i=1,d 
call putdec(STDOUT,a(i),10)
call putc(NEWLINE)
if (d<2) goto 10      # dimension too small
if (d>MPDIM) goto 20  # dimension too big
do i=2,d {
     if (a(i)<0) goto 30        # negative digit
     if (a(i)>=MPBASE) goto 40  # digit >= radix
     }
return
10 call putdec(STDOUT,code,5)
call putdec(STDOUT,d,10)
call putlin(s1,STDOUT)
goto 90
20 call putdec(STDOUT,code,5)
call putdec(STDOUT,d,10)
call putlin(s2,STDOUT)
goto 90
30 call putdec(STDOUT,code,5)
call putdec(STDOUT,d,10)
call putdec(STDOUT,i,10)
call putdec(STDOUT,a(i),10)
call putlin(s3,STDOUT)
goto 90
40 call putdec(STDOUT,code,5)
call putdec(STDOUT,d,10)
call putdec(STDOUT,i,10)
call putdec(STDOUT,a(i),10)
call putlin(s4,STDOUT)
90 call putc(NEWLINE)
stop
end

# MUL:  C = A*B      ssw  2/26/83
subroutine mul(a,b,c)
integer a(1), b(1), c(1), i,da,db,j,carry,x(2),k
# call mul1 if either a or b is single precision (dimension = 2)
if (a(1)<3) {call mul1(b,a(2),c); return}
if (b(1)<3) {call mul1(a,b(2),c); return}
MPTEST call mpchk(600,a)
MPTEST call mpchk(601,b)
da=a(1)
db=b(1)  # copy the dimensions into simple variables for do limits
do i=2,da
c(i)=0
do j=2,db {
     carry=0
     do i=2,da {
          k=i+j-2
          call muladd(a(i),b(j),c(k)+carry,x)
# x = a(i)*b(j) + c(i+j-2) + carry
          carry=x(2)
          c(k)=x(1)
          }
     k=da+j-1
     c(k)=carry
     }
c(1)=da+db-1
i=c(1)
if (c(i)==0) c(1) = c(1) - 1
MPTEST call mpchk(602,c)
return
end

# DIV:  Q = A/B; R = MOD(A,B)     ssw 4/4/83
subroutine div(a,b,q,r)
integer a(1),b(1),q(1),r(1),m,d,db,j,qhat,l,i
integer div1,y(2),x(2),borrow,tmp,k
#include 'mpcom.r' #define junk1,junk2,junk3 here as large enough scratch arrays
# file mpcom.r
common /mpcomm/ junkc(MPCDIM),junk1(MPDIM),junk2(MPDIM),junk3(MPDIM),junk4(MPDIM)
character junkc
integer junk1,junk2,junk3,junk4
CHAR*80 norm,q1,unn
norm = 'div norm failure'//char(0)
q1 = ' qhat - 1 '//char(0)
unn = 'div unnorm failure'//char(0)
# call div1 if  b  is single precision (dimension = 2)
if (b(1)<3) {r(1)=2; r(2)=div1(a,b(2),q); return}
MPTEST call mpchk(700,a)
MPTEST call mpchk(701,b) # check the input data
m=a(1)-b(1)+2  # dq will be m or m-1
db=b(1)
if (m>1) {  # the else is near the end of the subroutine.
     d=MPBASE/(1+b(db))  # step 1: normalize
     call mul1(a,d,junk1)  # u = junk1 = d*a
     call mul1(b,d,junk2)  # v = junk2 = d*b
     k=junk2(1); junk2(k+1)=0  # for use in step 6
     if (k!=db) {call putdec(STDOUT,db,10); call putdec(STDOUT,k,10) 
          call putc(NEWLINE); call error(norm)}
     l=junk1(1)  # l = da or da+1
     if (l==a(1)) {l=l+1; junk1(l)=0}
     for (j=0; j<m-1; j=j+1) {   # step 2: start j loop
          if (junk1(l)==junk2(db)) {   # step 3: find qhat
               qhat=MPBASE-1
               i=junk1(l-1)+junk2(db)}
          else {
#               t(1)=junk1(l-1)
#               t(2)=junk1(l)
#               i=divrem(t,junk2(db),qhat)} # qhat=t/junk2(db), i is dummy
		      call divrem(junk1(l),junk1(l-1),junk2(db),y)
		      i=y(2)	# maybe i is not a dummy after all
		      qhat=y(1)}
          call muladd(qhat,junk2(db-1),0,x)  # x = qhat*junk2(db-1) + 0
          while (x(2)>i|(x(2)==i&x(1)>junk1(l-2))){
               qhat=qhat-1
               MPTEST call putdec(STDOUT,qhat,6)
               MPTEST call putlin(q1,STDOUT)
               i=i+junk2(db)
               MPTEST call putdec(STDOUT,i,6)
               MPTEST call putc(NEWLINE)
               call muladd(qhat,junk2(db-1),0,x)
               }
          call mul1(junk2,qhat,junk3) # step 4: multiply and subtract
          if (junk2(1)==junk3(1)) {k=junk2(1)+1; junk3(1)=k; junk3(k)=0}
          borrow=0  # junk3 = qhat*junk2, perhaps padded with a high order zero
          k=l-junk2(1)+1
          for (i=2; i<=junk3(1); i=i+1) {  # subtract
               tmp=junk1(k)-junk3(i)-borrow
               if (tmp<0) {borrow=1; junk1(k)=tmp+MPBASE}  # borrow 1
               else       {borrow=0; junk1(k)=tmp}         # no borrow
               k=k+1
               }  # end i loop
          k=m-j         # step 5: test remainder
          if (borrow==0) {q(k)=qhat}
          else {q(k)=qhat-1
               borrow=0    # step 6: add back   'borrow' is 'carry' here
               k=l-junk2(1)+1
               for (i=2; i<=junk2(1); i=i+1) {
                    tmp=junk1(k)+junk2(i)+borrow
                    if (tmp<MPBASE) {borrow=0; junk1(k)=tmp}  # no carry
                    else {borrow=1; junk1(k)=tmp-MPBASE}      # carry 1
                    k=k+1
                    }  # ignore final carry
               } # end step 6
          l=l-1
          }  # step 7:  end loop on j
     while (db>2&junk1(db)==0) db=db-1  # condense junk1
     junk1(1)=db   # step 8: unnormalize
     i=div1(junk1,d,r)  # i had better be zero
     if (i!=0) {call putdec(STDOUT,i,10); call putc(NEWLINE); call error(unn)}
     if (q(m)==0 & m>2) {q(1)=m-1} else {q(1)=m}  # set dimension of q
     }  else  {   # else m <= 1   the 'if' was near beginning of subroutine
     call mpcopy(a,r)  # set r=a
     q(1)=2  # q=0
     q(2)=0
     } 
MPTEST call mpchk(702,q)
MPTEST call mpchk(703,r)
return
end
# file: mp2.r;  author:  ssw;  modified by jwt so output can be to any file
# MP subroutines MPITOC, MPCTOI, MPTDEC, MPITOH, MPHTOI, MPTHEX,
#                MPCMP, MPSQRT, MPODIV, MPGOCD, MPOWER, JACOBI
# UPPER CASE subroutines are for general use.
# lower case subroutines are mostly for private internal use only.

# MPITOC:  length = MPITOC(mpint,string,size);  convert mp integer  mpint
#   to character string  string  of decimal digits of length MPITOC <= size.
#   ssw 3/29/83
integer function mpitoc(int,str,size)
integer int(1),size,div1,d,i,j
character str(1),digits(10),k
CHAR*80 big
#include 'mpcom.r'  # define junk1 here as large enough scratch array
# file mpcom.r
common /mpcomm/ junkc(MPCDIM),junk1(MPDIM),junk2(MPDIM),junk3(MPDIM),junk4(MPDIM)
character junkc
integer junk1,junk2,junk3,junk4
data digits/'0','1','2','3','4','5','6','7','8','9'/
big = 'int too big for mpitoc'//char(0)
if (int(1) > MPDIM) call error(big)
call mpcopy(int,junk1)
str(1)=EOS
i=1
repeat {                # generate digits
     i=i+1
     d=div1(junk1,10,junk1) # d=mod(junk1,10); junk1=junk1/10
     str(i)=digits(d+1)
     } until ((junk1(1)==2&junk1(2)==0)|i>=size) # first cond is 'junk1=0'
mpitoc=i-1
for (j=1; j<i; j=j+1) {     # then reverse
     k=str(i)
     str(i)=str(j)
     str(j)=k
     i=i-1
     }
return
end

# MPCTOI:  convert string of decimal digits at  in(i)  to mp integer  out  and 
#          increment  i.    ssw 3/29/83
subroutine mpctoi(in,i,out)
character in(1)
integer index, d, i, out(1), x(2)
equivalence (x(2),d)
CHAR*10 digits
data digits /'0123456789'/
data x(1)/2/
while (in(i)==BLANK)      
i=i+1
out(1)=2
for (out(2)=0; in(i)!=EOS; i=i+1) {
     d=index(digits,in(i))
     if (d==0)   # non-digit
     break
     call mul1(out,10,out)    # out=10*out
     d=d-1
     call add(out,x,out)      # out=out+d-1
     }
return
end

# MPTDEC:  put mp integer  n  as decimal digits into field of width >= w.
#  ssw 3/29/83
subroutine mptdec(lunit,n,w)
integer mpitoc,i,n(1),nd,w,lunit,j,fputc
#include 'mpcom.r'   # define junkc in this common block
# file mpcom.r
common /mpcomm/ junkc(MPCDIM),junk1(MPDIM),junk2(MPDIM),junk3(MPDIM),junk4(MPDIM)
character junkc
integer junk1,junk2,junk3,junk4
nd=mpitoc(n,junkc,MPCDIM)
for (i=nd+1; i<=w; i=i+1)
j= fputc(lunit,BLANK)
for (i=1; i<=nd; i=i+1)
j= fputc(lunit,junkc(i))
return
end

# MPCMP: value = +1 if A>B, = 0 if A=B, and = -1 if A<B, i.e., value = sgn(A-B)
#  ssw 5/21/83
integer function mpcmp(a,b)
integer a(1), b(1), i
if (a(1)<b(1)) {mpcmp=-1}  # compare dimensions first
else if (a(1)>b(1)) {mpcmp=1}
else {  # a(1) = b(1), so compare the digits from high to low order
     mpcmp=0  # assume a = b until proven otherwise
     for (i=a(1); i>1; i=i-1) {
          if (a(i)<b(i)) {mpcmp=-1; break}
          else if (a(i)>b(i)) {mpcmp=1; break}
          }
     }
return
end

# MPODIV: Remove largest power (2**L) of 2 from A (so A is both I/O).
# Function value MPODIV = L.  But MPODIV = -1 if A = 0.  ssw 6/6/83
integer function mpodiv(a)
integer a(1),div1,i,d,j,l
CHAR*80 e
e = 'error in MPODIV'//char(0)
if (a(1)==2&a(2)==0) {mpodiv=-1; return}
# first remove the largest power of MPBASE
for (i=2; a(i)==0; i=i+1) ;  # loop has null scope; it just the first a(i)!=0
i=i-2
if (i>0) {d=a(1)-i
     do j=2,d {
          l=j+i
          a(j)=a(l)
          }
     a(1)=d
     }
# now remove a power of two smaller than MPBASE
j=0
d=a(2)
while (mod(d,2)==0) {j=j+1; d=d/2}
mpodiv=MPBITS*i+j
if (j>0) {
     if (div1(a,2**j,a)!=0) { # there is an error if the remainder is non-zero
          call mpchk(550,a)
          call error(e)
          }
     }
return
end

# MPGOCD: C = greatest odd common divisor of A and B   ssw 6/6/83
subroutine mpgocd(a,b,c)
integer a(1),b(1),c(1),i,j,mpodiv
#include 'mpcom.r'  # defines junk1 as a large scratch array
# file mpcom.r
common /mpcomm/ junkc(MPCDIM),junk1(MPDIM),junk2(MPDIM),junk3(MPDIM),junk4(MPDIM)
character junkc
integer junk1,junk2,junk3,junk4
call mpcopy(a,c)
call mpcopy(b,junk1)
i=mpodiv(c)  # remove all powers of 2 from c
if (mpodiv(junk1)<0) return  # check for zero b
if (i<0) {call mpcopy(junk1,c); return} # check for zero a.
i=mpcmp(c,junk1)  # compare c with junk1
while (i!=0) { # repeat until c = junk1
     if (i<0) {  # if c < junk1
          call sub(junk1,c,junk1)  # junk1 = junk1 - c
          j=mpodiv(junk1)   # remove 2's
          }
     else {  # else c > junk1
          call sub(c,junk1,c)  # c = c - junk1
          j=mpodiv(c)    # remove 2's
          }
     i=mpcmp(c,junk1)
     }
return
end

# MPOWER: Compute A = B**E (mod M)  Assume 0 < B < M.  ssw 6/9/83
subroutine mpower(b,e,m,a)
integer b(1),e(1),m(1),a(1),bits(MPBITS),i,n,ebit
#include 'mpcom.r'  # define junkc and junk4 as large scratch arrays.
 # junk1, junk2, junk3 are not used because div uses them and we call div.
# file mpcom.r
common /mpcomm/ junkc(MPCDIM),junk1(MPDIM),junk2(MPDIM),junk3(MPDIM),junk4(MPDIM)
character junkc
integer junk1,junk2,junk3,junk4
call mpcopy(b,a)  # initialize a = b
for (n=e(1); n>1; n=n-1) {  # begin outer loop
     ebit=e(n)
     i=0
     do j=1,MPBITS {  # break ebit = e(n) into bits, saved in bits
          i=i+1
          bits(i)=mod(ebit,2)   # or use and(ebit,1) if available
          ebit=ebit/2           # or shift ebit one bit
          if (n==e(1)&ebit==0) {i=i-1; break}
          }
     for (; i>0; i=i-1) {  # i = MPBITS except first time
          call mul(a,a,junkc)  # junkc = a*a
          call div(junkc,m,junk4,a)  # a = junkc mod m; junk4 is ignored
          if (bits(i)!=0) {
               call mul(a,b,junkc)  # junkc = a*b
               call div(junkc,m,junk4,a)  # a = junkc mod m; ignore junk4
               }
          }
     }
return
end

# JACOBI: Compute Jacobi symbol (pa/qa)   ssw 6/6/83
# We had better have qa odd, qa > 0, and pa > 0.
integer function jacobi(pa,qa)
integer pa,qa,p,q,r,itmp
CHAR*80 e
e = 'bad argument for JACOBI'//char(0)
jacobi=1
if (qa<=0|pa<0|mod(qa,2)==0) {  # we check the arguments
     call putdec(STDOUT,pa,5)
     call putdec(STDOUT,qa,5)
     call putc(NEWLINE)
     call error(e)
     }
if (qa==1|pa==1) return  # value is 1 if either pa or qa = 1
p=pa  # copy over pa and qa so we can change them
q=qa  # without changing the originals
if (q<=p) p=mod(p,q)
if (p==0) {jacobi=0; p=1}
while (p!=1) {  # repeat until p=1
	itmp=1+and(p,3)
#     switch (mod(p,4))  {  # break into 4 cases of p mod 4
	switch (itmp) {
          case 1: p=p/4
          case 2: {r=mod(q,p)
               if (r==0) {p=1; jacobi=0} # exit if r = 0
               else {q=p; p=r}
               }
          case 3: {p=p/2
               r=mod(q,8)  # or use r=and(q,7) if available
               if (r==3|r==5) jacobi=-jacobi
               }
          case 4: {if (mod(q,4)==3) jacobi=-jacobi
               r=mod(q,p)
               if (r==0) {p=1; jacobi=0}  # exit if r = 0
               else {q=p; p=r}
               }
          }
     }
return
end
